/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package distribution;

import cern.jet.random.AbstractDiscreteDistribution.*;
import cern.jet.random.Empirical;
import cern.jet.random.EmpiricalWalker;
import cern.jet.random.engine.RandomEngine;
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.logging.Level;
import java.util.logging.Logger;
import umontreal.iro.lecuyer.probdist.EmpiricalDist;
import pedigreereconstruction.*;

//import org.apache.commons.math3.random.EmpiricalDistribution;
// another possiblility: umontreal.iro.lecuyer.probdist.DiscreteDistribution
/**
 *
 * @author arthiramachandran
 */
public class RelationshipDistribution {

    public static final int COLT = 1, UMONTREAL = 2, APACHECOMMONS = 3;
    public static final int statsPackage = UMONTREAL;
    String relationshipName;
    //public EmpiricalDistribution lengthDistribution = null;
    //public EmpiricalDistribution countDistribution = null;
    // UMONTREAL
    public EmpiricalDist lengthDistribution = null;
    public EmpiricalDist countDistribution = null;
    //COLT
    //public EmpiricalWalker lengthDistribution;
    //public EmpiricalWalker countDistribution;
    public ArrayList<Double> lengths = new ArrayList<Double>();
    public ArrayList<Double> counts = new ArrayList<Double>();
    Double maxActualLength = Double.NEGATIVE_INFINITY,
            maxActualCount = Double.NEGATIVE_INFINITY;
    // some sort of distribution of means, parameters
    public boolean symmetric;

    // TODO: figure out a way to compute cdf with doubles rather than int
    // TODO: is this correct to handle empirical estimates which don't have enough info?
    /*public double calculateCDF(Double length, int count) {
        double normalization = 0;
        for (RelationshipDistribution d : PedigreeReconstruction.knownPairwiseDistributions.values()) {
            normalization += d.calculateCDFSub(length, count);
        }

        return this.calculateCDFSub(length, count) / normalization;
    }
     * 
     */

    public double calculateCDF(Double length, int count) {
        //if (this.relationshipName.equals("GRANDPARENT_GRANDCHILD"))
        //    return 0;
        //System.out.println("cdfcalc");
        if (lengths.size() > 1) {
            //System.out.println(relationshipName + " " + length + " " + count);
            /*
            System.out.println("&&& " + lengthDistribution.getNumericalVariance() + "  "
            + lengthDistribution.getBinCount() + " "
            + lengthDistribution.cumulativeProbability(length) + "  "
            + countDistribution.getBinCount() + " "
            + countDistribution.cumulativeProbability(count));
             * 
             */
            /*
            System.out.println(length + " "
            + lengthDistribution.cdf(length) + " "
            + lengthDistribution.getSampleMean() + " "
            + lengthDistribution.getSampleStandardDeviation() + " "
            + count + " " + countDistribution.cdf(count) + " "
            + countDistribution.getSampleMean() + " "
            + countDistribution.getSampleStandardDeviation());
             * 
             */


            //System.out.println(lengthDistribution.getSampleStats());
            //System.out.println(countDistribution.getSampleStats());
            //System.out.println("--" + lengthDistribution.cumulativeProbability(length));
            //System.out.println(countDistribution.cumulativeProbability(count));
            if (statsPackage == UMONTREAL) {
                double lengthPdf =
                        lengthDistribution.cdf(length + 200) - lengthDistribution.cdf(length);
                double countPdf =
                        countDistribution.cdf(count + 5) - countDistribution.cdf(count);
                /*System.out.println("$$$$ " + this.relationshipName + " " + length + " " + maxActualLength + " "
                + lengthDistribution.cdf(length) + " " + lengthDistribution.cdf(maxActualLength)
                + " " + lengthDistribution.cdf(length + 100) + " " + lengthPdf);
                System.out.println("$$$$ " + this.relationshipName + " " + count + " " + maxActualCount + " "
                + countDistribution.cdf(count) + " " + countDistribution.cdf(maxActualCount)
                + " " + countDistribution.cdf(count + 5) + " " + countPdf);
                System.out.println(lengthPdf * countPdf);
                 * 
                 */

                // HACK to account for really high sharing pairs
                if (length > maxActualLength && this.relationshipName.equals("PARENT_CHILD")) {
                    lengthPdf = 1;
                } else if (length > maxActualLength) {
                    lengthPdf =
                            (lengthDistribution.cdf(length) - lengthDistribution.cdf(maxActualLength)) / (length - maxActualLength) * 100;
                    /*System.out.println("$$$$ " + this.relationshipName + " " + length + " " + maxActualLength + " "
                    + lengthDistribution.cdf(length) + " " + lengthDistribution.cdf(maxActualLength));
                     * 
                     */
                }
                if (count > maxActualCount && this.relationshipName.equals("PARENT_CHILD")) {
                    countPdf = 1;
                }

                int i = 200;
                while (lengthPdf == 0) {
                    lengthPdf =
                            (lengthDistribution.cdf(length + i) - lengthDistribution.cdf(length))
                            / i * 100;
                    i += 100;
                    /*System.out.println("$$$$## " + length + " " + i + "-" + (length + i) + " "
                    + lengthDistribution.cdf(length) + " " + lengthDistribution.cdf(maxActualLength)
                    + " " + lengthDistribution.cdf(length + i));
                     * 
                     */
                }
                if (countPdf == 0) {
                    countPdf =
                            lengthDistribution.cdf(count) - lengthDistribution.cdf(maxActualCount);
                    System.out.println("$$$$## " + count + " " + maxActualCount);
                }





                //System.out.println(lengthPdf * countPdf);
                if (lengthPdf * countPdf == 0) {
                    System.out.println("AHH!!");
                    System.out.println("$$$$ " + this.relationshipName + " " + length + " " + maxActualLength + " "
                            + lengthDistribution.cdf(length) + " " + lengthDistribution.cdf(maxActualLength)
                            + " " + lengthDistribution.cdf(length + 200) + " " + lengthPdf);
                    System.out.println("$$$$ " + this.relationshipName + " " + count + " " + maxActualCount + " "
                            + countDistribution.cdf(count) + " " + countDistribution.cdf(maxActualCount)
                            + " " + countDistribution.cdf(count + 5) + " " + countPdf);
                    System.out.println(lengthPdf * countPdf);
                }


                return (lengthPdf * countPdf);
            } else if (statsPackage == COLT) {
                System.out.println("Wrong package!!");
                return (lengthDistribution.cdf(length.intValue()) * countDistribution.cdf(count));

            } else {
                System.out.println("Temp fix");
                return 1;
            }
            //return (lengthDistribution.cumulativeProbability(length) * countDistribution.cumulativeProbability(count));
        } else {
            System.out.println("Relationship " + relationshipName + " doesn't have enough information to make correct calculations");
            return 1;
        }
    }

    public RelationshipDistribution(String name, boolean symmetric,
            ArrayList<Double> lengthPdf,
            ArrayList<Double> countPdf) throws Exception {
        this.relationshipName = name;
        this.symmetric = symmetric;

        maxActualLength = getMaxValue(lengthPdf);
        maxActualCount = getMaxValue(countPdf);


        //addLaplacianCorrection(lengthPdf, countPdf);
        // Laplacian correction
        for (int i = 1; i < 100; i++) {
            lengthPdf.add(60.0 * i); // 6000cM possible sharing
            countPdf.add(1.5 * i); // todo: random number of max total shared segments
        }

        /*
        if (statsPackage == COLT) {
        this.lengthDistribution = new EmpiricalWalker(lengthPdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        this.countDistribution = new EmpiricalWalker(countPdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        }
         * 
         */

        /*
        lengthDistribution = new EmpiricalDistribution(5);
        lengthDistribution.load(lengthPdf);
        countDistribution = new EmpiricalDistribution(5);
        countDistribution.load(lengthPdf);
         * 
         */
        for (int i = 0; i < lengthPdf.size(); i++) {
            lengths.add(lengthPdf.get(i));
            counts.add(countPdf.get(i));
        }
        int nBins = lengthPdf.size();
        double[] lpdf = new double[nBins];
        double[] cpdf = new double[nBins];

        for (int i = 0; i < lengthPdf.size(); i++) {
            lpdf[i] = lengthPdf.get(i);
            cpdf[i] = countPdf.get(i);
        }

        // sort the arrays of data
        if (statsPackage == UMONTREAL) {
            Arrays.sort(lpdf);
            Arrays.sort(cpdf);
            lengthDistribution = new EmpiricalDist(lpdf);
            countDistribution = new EmpiricalDist(cpdf);
        }
        // todo: this is stupid
    }

    public RelationshipDistribution(String name, boolean symmetric,
            String lengthCountPdfFileName)
            throws Exception {
        this.relationshipName = name;
        this.symmetric = symmetric;

        System.out.println("Reading " + lengthCountPdfFileName);

        ArrayList<Double> lengthPdf = new ArrayList<Double>();
        ArrayList<Double> countPdf = new ArrayList<Double>();
        BufferedReader in = null;
        try {
            in = new BufferedReader(new FileReader(lengthCountPdfFileName));
            String line;
            while ((line = in.readLine()) != null) {
                if (line.startsWith("id")) // header line
                {
                    continue;
                }
                String[] fields = line.split(",");
                try {
                    lengthPdf.add(Double.parseDouble(fields[3]));
                    countPdf.add(Double.parseDouble(fields[2]));
                } catch (NumberFormatException e) {
                    System.out.println(
                            "Skipping malformed line in file "
                            + lengthCountPdfFileName + "\n Malformed line:\n"
                            + line
                            + "\nAll lengths and counts of IBD need to be doubles");
                    continue;
                }
            }
        } catch (IOException ex) {
            Logger.getLogger(RelationshipDistribution.class.getName()).log(
                    Level.SEVERE, null, ex);
        } finally {
            try {
                in.close();
            } catch (IOException ex) {
                Logger.getLogger(RelationshipDistribution.class.getName()).log(
                        Level.SEVERE, null, ex);
            }
        }

        if (lengthPdf.size() != countPdf.size()) {
            throw new Exception(
                    "Input file " + lengthCountPdfFileName + " has "
                    + "incomplete information; lengths and counts of shared "
                    + "IBD are inconsistent (the number of individuals with "
                    + "valid counts is different than the number of "
                    + "individuals with valid lengths");
        }
        //System.out.println("Length: " + lengthPdf);
        //System.out.println("Counts: " + countPdf);
        maxActualLength = getMaxValue(lengthPdf);
        maxActualCount = getMaxValue(countPdf);

        // Laplacian correction - is it??
        //addLaplacianCorrection(lengthPdf, countPdf);
        for (int i = 1; i < 100; i++) {
            lengthPdf.add(60.0 * i); // 6000cM possible sharing
            countPdf.add(1.5 * i); // todo: random number of max total shared segments
        }

        int nBins = lengthPdf.size();
        double[] lpdf = new double[nBins];
        double[] cpdf = new double[nBins];
        /* need to 
         *  sort
         *  put into bins
         * when computing pdf
         *  figure out which bin
         *  and get the pdf/cdf
         */


        for (int i = 0; i < lengthPdf.size(); i++) {
            lpdf[i] = lengthPdf.get(i);
            //System.out.println(lpdf[i]);
            cpdf[i] = countPdf.get(i);
        }

        //double[] lengthPdf  = new double[]();
        /*
        if (statsPackage == COLT) {
        this.lengthDistribution = new EmpiricalWalker(lpdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        this.countDistribution = new EmpiricalWalker(cpdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        }
         * 
         */
        // TODO: decent metric?
        /*
        lengthDistribution = new EmpiricalDistribution(lengthPdf.size() / 10);
        lengthDistribution.load(lpdf);
        countDistribution = new EmpiricalDistribution(countPdf.size() / 10);
        countDistribution.load(cpdf);
         * 
         */
        if (statsPackage == UMONTREAL) {
            Arrays.sort(lpdf);
            Arrays.sort(cpdf);

            lengthDistribution = new EmpiricalDist(lpdf);
            countDistribution = new EmpiricalDist(cpdf);
        }
        lengths = lengthPdf;
        counts = countPdf;
    }
    /* 
     * make an abstract Distribution class
     */

    public String getName() {
        return this.relationshipName;
    }

    public double getRandomRelationshipProbability() {
        // TODO: better random?
        /*System.out.println("Randoms: " + random1 + " " + random2 + " " +lengthDistribution.cdf(lengthDistribution.getObs(random1))
         * countDistribution.cdf(countDistribution.getObs(random2)));
         * 
         */
        if (statsPackage == UMONTREAL) {
            int random1 =
                    (int) Math.floor((Math.random() * lengthDistribution.getN()));
            int random2 =
                    (int) Math.floor((Math.random() * countDistribution.getN()));
            return lengthDistribution.cdf(lengthDistribution.getObs(random1))
                    * countDistribution.cdf(countDistribution.getObs(random2));
        }
        /*
        if (statsPackage == COLT) {
        return lengthDistribution.cdf(lengthDistribution.nextInt())
         * countDistribution.cdf(countDistribution.nextInt());
        }
         * 
         */
        return 1;
        /*return lengthDistribution.cumulativeProbability(lengthDistribution.sample())
         * countDistribution.cumulativeProbability(countDistribution.sample());
         * 
         */
    }

    public void addDataPoint(int numberOfSharedSegments, double totalSharedIBD) {
        this.lengths.add(totalSharedIBD);
        this.counts.add((double) numberOfSharedSegments);
        double[] lpdf = new double[lengths.size()];
        double[] cpdf = new double[counts.size()];
        for (int i = 0; i < lengths.size(); i++) {
            lpdf[i] = lengths.get(i);
            cpdf[i] = counts.get(i);
        }
        // TODO: check this!!!
        if (statsPackage == UMONTREAL) {
            Arrays.sort(lpdf);
            Arrays.sort(cpdf);

            this.lengthDistribution = new EmpiricalDist(lpdf);
            this.countDistribution = new EmpiricalDist(cpdf);
        }
        /*
        this.lengthDistribution.load(lpdf);
        this.countDistribution.load(cpdf);
         * 
         */
        /*
        if (statsPackage == COLT) {
        this.lengthDistribution = new EmpiricalWalker(lpdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        this.countDistribution = new EmpiricalWalker(cpdf,
        Empirical.NO_INTERPOLATION, RandomEngine.makeDefault());
        }
         * 
         */
    }

    public int numberOfDataPoints() {
        return this.lengths.size();
    }

    private void addLaplacianCorrection(double[] lengthPdf, double[] countPdf) {
        throw new UnsupportedOperationException("Not yet implemented");
    }

    public String getMinValues() {
        double cmin = Double.POSITIVE_INFINITY, lmin = Double.POSITIVE_INFINITY;

        for (double c : counts) {
            if (c < cmin) {
                cmin = c;
            }
        }
        for (double l : lengths) {
            if (l < lmin) {
                lmin = l;
            }
        }
        return (cmin + "," + lmin);
    }

    public String getMaxValues() {
        double cmax = Double.NEGATIVE_INFINITY, lmax = Double.NEGATIVE_INFINITY;

        for (double c : counts) {
            if (c > cmax) {
                cmax = c;
            }
        }
        for (double l : lengths) {
            if (l > lmax) {
                lmax = l;
            }
        }
        return (maxActualLength + ", " + maxActualCount);
        //return (cmax + "," + lmax);
    }

    public String detailedCalculateCdf(double length,
            int count) {

        System.out.println("Detail");
        double lengthPdf =
                lengthDistribution.cdf(length + 100) - lengthDistribution.cdf(length);
        double countPdf =
                countDistribution.cdf(count + 5) - countDistribution.cdf(count);
        System.out.println("$$$$ " + this.relationshipName + " " + length + " " + maxActualLength + " "
                + lengthDistribution.cdf(length) + " " + lengthDistribution.cdf(maxActualLength)
                + " " + lengthDistribution.cdf(length + 100) + " " + lengthPdf);
        System.out.println("$$$$ " + this.relationshipName + " " + count + " " + maxActualCount + " "
                + countDistribution.cdf(count) + " " + countDistribution.cdf(maxActualCount)
                + " " + countDistribution.cdf(count + 5) + " " + countPdf);
        System.out.println(lengthPdf * countPdf);

        //System.out.println(this.getMaxValues());
        return (length + "=" + lengthPdf + ", " + count + "=" + countPdf);
    }

    private Double getMaxValue(ArrayList<Double> list) {
        double max = Double.NEGATIVE_INFINITY;

        for (double element : list) {
            if (element > max) {
                max = element;
            }
        }
        return max;
    }
}
